/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.materialMedia.constitutiveModel;

import java.io.Serializable;
import java.util.HashMap;

import br.ufmg.dees.insane.analysisModel.AnalysisModel;
import br.ufmg.dees.insane.analysisModel.PlaneStrain;
import br.ufmg.dees.insane.analysisModel.PlaneStress;
import br.ufmg.dees.insane.analysisModel.PolarPlaneStrain;
import br.ufmg.dees.insane.analysisModel.PolarPlaneStress;
import br.ufmg.dees.insane.continuousPointModel.ContinuousPointModel;
import br.ufmg.dees.insane.materialMedia.material.Material;
import br.ufmg.dees.insane.materialMedia.material.VonMises;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;
import br.ufmg.dees.insane.util.LinearEquationSystems;

/** A class representing ElastoPlasticConstModel
 * @author Jamile Fuina
 * @since 22/05/2005
 */
public abstract class ElastoPlasticConstModel implements ConstitutiveModel ,Serializable {
	
	private static final long serialVersionUID = 1L;

    /** The tolerance used to verify plastification.*/ 
    private double EP_TOL = 1.0e-10; 
    
    /** The number maximum of iteractions.*/ 
    private int EP_MAX_ITE = 50; 
    
    /** The plastic multiplier used in associative plasticity.*/
    public final String PLASTIC_MULTIPLIER = "PLASTIC_MULTIPLIER";
    
    /** The algorithmic moduli used in associative plasticity.*/
    public final String ALGORIT_MOD = "ALGORIT_MOD";
    
    /** The derivative of yield surface w.r.t. stress used in associative plasticity.*/
    public final String dFdSIG = "dFdSIG";
    
    /** The derivative of yield surface w.r.t. internal variables used in associative plasticity.*/
    public final String dFdq = "dFdq";
    
    /** The number {dF/dSIG dF/dQ}*[A]*{dF/dSIG dF/dQ} used in associative plasticity.*/
    public final String DENOM = "DENOM";
    
    /** The plastic deformation used in associative plasticity.*/
    public final String PLASTIC_DEFORM = "PLASTIC_DEFORM";
    
    /** The iteractive plastic deformation used in associative plasticity.*/
    public final String ITE_PLASTIC_DEFORM = "ITE_PLASTIC_DEFORM";
    
    /** The dimension of Hessian of Lagrangian function.*/
    protected int dimA ;
    
    /**Constructs this Constitutive Model.*/
    public ElastoPlasticConstModel()
	{
	}
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountCt(br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IMatrix mountCt(AnalysisModel am, Material mat, HashMap hm) {
    	this.check(am,mat);
        if (Math.abs(((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue()) <= EP_TOL){
            return (am.getDualInternalVariableOperator(mat.getPt()));
        }
        IMatrix Ct = new IMatrix(am.getNumberOfDualInternalVariables(),am.getNumberOfDualInternalVariables());
        IVector aux = new IVector(am.getNumberOfDualInternalVariables());
        aux.mul((IMatrix)hm.get(this.ALGORIT_MOD), (IVector)hm.get(this.dFdSIG));
        for(int i=0; i<am.getNumberOfDualInternalVariables(); i++)
            for(int j=0; j<am.getNumberOfDualInternalVariables(); j++)
                Ct.setElement(i,j,((IMatrix)hm.get(this.ALGORIT_MOD)).getElement(i,j) - (aux.getElement(i) * aux.getElement(j)/((Double)hm.get(this.DENOM)).doubleValue()));         
        return Ct;
     }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountCs(br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IMatrix mountCs(AnalysisModel am, Material mat, HashMap hm) {
        this.check(am,mat);
        return (am.getDualInternalVariableOperator(mat.getPs()));
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountDualInternalVariableVector(br.ufmg.dees.insane.util.IVector, br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IVector mountDualInternalVariableVector(IVector e, AnalysisModel am, Material mat, HashMap<Object, Object> mp) {
    	this.check(am,mat);
        IMatrix aux1 = am.getDualInternalVariableOperator(mat.getPs());
        
        IMatrix aux2 = aux1.inverse();
        IVector stress= new IVector(am.getNumberOfDualInternalVariables());
        try{
            this.ClosestPoint(e,aux1,aux2, am, mat, mp, stress);
        }
        catch(ReturnNotAchieved exc){
            System.out.println(exc.getMessage() + "to try with cutting-plane");
            try{
                this.CuttingPlane(e,aux1, am, mat, mp, stress);  
            }
            catch(ReturnNotAchieved ex){
                System.out.println(ex.getMessage() + "the system exit"); 
                System.exit(1);
        }   
            }
            return stress;
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountDualInternalVariableVector(br.ufmg.dees.insane.util.IVector, br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IVector mountDualInternalVariableVector(IVector e, ContinuousPointModel am, Material mat, HashMap mp) {
    	return null;
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountC(br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IMatrix mountC(AnalysisModel am, Material mat, HashMap mp) {
    	this.check(am,mat);
        return (am.getDualInternalVariableOperator(mat.getPs()));        
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#mountC(br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material, java.util.HashMap)
     */
    public IMatrix mountC(ContinuousPointModel am, Material mat, HashMap mp) {
    	return null;
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#check(br.ufmg.dees.insane.analysisModel.AnalysisModel)
     */
    public void check(AnalysisModel am) throws AnalysisModelNotValid {
        if (!(am instanceof PlaneStrain ||
            am instanceof PlaneStress ||
            am instanceof PolarPlaneStrain ||
            am instanceof PolarPlaneStress)) {
            throw new AnalysisModelNotValid(this.getLabel()+" and "+am.getLabel());
        }
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#check(br.ufmg.dees.insane.analysisModel.AnalysisModel)
     */
    public void check(ContinuousPointModel am) throws ContinuousPointModelNotValid {
    }
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#check(br.ufmg.dees.insane.materialMedia.material.Material)
     */
    public void check(Material mat) throws MaterialNotValid {
        if (!(mat instanceof VonMises)) {
            throw new MaterialNotValid(this.getLabel()+" and "+mat.getClass());
        }
    }
 
    /**Returns the state of stress from a given state of strain using the Closest-Point Projection Algorithm. Evaluates the plastic strains and the internal variables. 
     * @param et The total strain.
     * @param c  The elastic constitutive matrix. 
     * @param c1  The inverse of elastic constitutive matrix.
     * @return Returns the stress vector.
     */
    private boolean ClosestPoint (IVector et, IMatrix c, IMatrix c1, AnalysisModel am, Material mat, HashMap<Object, Object> hm, IVector stress) throws ReturnNotAchieved
    {
        boolean conv = true;
        IVector ee = new IVector(am.getNumberOfDualInternalVariables());
        IVector r = new IVector(this.dimA);
        IVector v = new IVector(this.dimA);
        IVector ep0 = ((IVector)hm.get(this.ITE_PLASTIC_DEFORM));
        IVector alpha0 = this.getInternalVar(mat, hm);
        IVector ep = ep0;
        IVector alpha = alpha0;
        LinearEquationSystems sel = new LinearEquationSystems();
        int k = 0;
        do{
          ee.sub(et,ep);
          stress.mul(c,ee);
          IVector q = internalVarStress(alpha, mat);
          double f = this.yieldFunction(stress,q,am,mat);
          hm.put(this.dFdSIG,yieldGradStress (stress,q,am,hm));
          hm.put(this.dFdq ,yieldGradIntVar (stress,q,mat,hm));
          for(int i=0; i<am.getNumberOfDualInternalVariables(); i++)               
              r.setElement(i,-ep.getElement(i) + ep0.getElement(i) + ((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue()*((IVector)hm.get(this.dFdSIG)).getElement(i));          
          for(int i=0; i<this.getNumInternalVar(mat); i++)             
              r.setElement(i + am.getNumberOfDualInternalVariables(),-alpha.getElement(i) + alpha0.getElement(i) + ((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue()*((IVector)hm.get(this.dFdq)).getElement(i));          
          double err = r.norm();
          if((f <= EP_TOL) && (err <= EP_TOL)){
        	  break; 	  
          }
          if(k > EP_MAX_ITE){
              conv = false;
              throw new ReturnNotAchieved("closest-point");
          }
          IMatrix D = plasticModuli(alpha, mat);
          IMatrix D1 = D.inverse();
          IMatrix A = yieldHessian(stress,q,am);
          A.scale(((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue());      
          for(int i=0; i<am.getNumberOfDualInternalVariables(); i++){
              for(int j=0; j<am.getNumberOfDualInternalVariables(); j++){
                 A.setElement(i,j, A.getElement(i,j) + c1.getElement(i,j));
                 hm.put(this.ALGORIT_MOD,A);
              }
          }
          for(int i=0; i<this.getNumInternalVar(mat); i++)
              for(int j=0; j<this.getNumInternalVar(mat); j++)
                 A.setElement(i + am.getNumberOfDualInternalVariables(),j + am.getNumberOfDualInternalVariables(), A.getElement(i + am.getNumberOfDualInternalVariables(),j + am.getNumberOfDualInternalVariables()) + D1.getElement(i,j));
          for(int i=0; i<am.getNumberOfDualInternalVariables(); i++) 
              v.setElement(i,((IVector)hm.get(this.dFdSIG)).getElement(i));
          for(int i=0; i<this.getNumInternalVar(mat); i++) 
              v.setElement(i + am.getNumberOfDualInternalVariables(),((IVector)hm.get(this.dFdq)).getElement(i));
          IVector aux1 = sel.solverX(A,v);
          Double denom = new Double(v.dot(aux1));
          hm.put(this.DENOM,denom);
          IVector aux2 = sel.solverX(A,r);
          double ddg = (f - aux2.dot(v))/((Double)hm.get(this.DENOM)).doubleValue();
          for(int i=0; i<am.getNumberOfDualInternalVariables(); i++) 
              v.setElement(i, r.getElement(i) + ddg * ((IVector)hm.get(this.dFdSIG)).getElement(i));
          for(int i=0; i<this.getNumInternalVar(mat); i++) 
              v.setElement(i + am.getNumberOfDualInternalVariables(),r.getElement(i) + ddg * ((IVector)hm.get(this.dFdq)).getElement(i));
          v = sel.solverX(A,v);          
          double dep;        
          for(int i=0; i<am.getNumberOfDualInternalVariables(); i++){
              dep = 0.0;
              for(int j=0; j<am.getNumberOfDualInternalVariables(); j++)
                  dep += c1.getElement(i,j)* v.getElement(j);
              ep.setElement(i,ep.getElement(i) + dep);
          }          
          double dalpha;        
          for(int i=0; i<this.getNumInternalVar(mat); i++){
              dalpha = 0.0;
              for(int j=0; j<this.getNumInternalVar(mat); j++)
                  dalpha += D1.getElement(i,j)* v.getElement(j + am.getNumberOfDualInternalVariables());
              alpha.setElement(i,alpha.getElement(i) + dalpha);
          }          
          Double inc = new Double(((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue()+ ddg);
          hm.put(this.PLASTIC_MULTIPLIER,inc);           
          k++;
        }
        while(true);
        if(conv && ((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue() > 0.0){
            IMatrix aux3 =((IMatrix)hm.get(this.ALGORIT_MOD)).inverse();
            hm.put(this.ALGORIT_MOD, aux3);             
            this.setPlasticVar(ep, alpha, mat, hm);
        }
        else hm.put(this.ALGORIT_MOD,c);
        return conv;
    }
    
    /**Returns the state of stress from a given state of strain using the CuttingPlane Plane Algorithm. Evaluates the plastic strains and the internal variables. 
     * @param et The total strain.
     * @param c  The elastic constitutive matrix. 
     * @param c1  The inverse of elastic constitutive matrix.
     * @return Returns the stress vector.
     */
    private boolean CuttingPlane (IVector et, IMatrix c, AnalysisModel am, Material mat, HashMap<Object, Object> hm, IVector stress) throws ReturnNotAchieved
    {
        boolean conv = true;
        this.dimA = am.getNumberOfDualInternalVariables() + this.getNumInternalVar(mat);        
        IVector ee = new IVector(am.getNumberOfDualInternalVariables());        
        IVector w = new IVector(this.dimA);
        IVector v = new IVector(am.getNumberOfDualInternalVariables());
        IVector ep0 = ((IVector)hm.get(this.ITE_PLASTIC_DEFORM));
        IVector alpha0 = this.getInternalVar(mat, hm);
        IVector ep = ep0;
        IVector alpha = alpha0;        
        int k = 0;
        do{
          ee.sub(et,ep);
          stress.mul(c,ee);
          IVector q = internalVarStress(alpha,mat);
          double f = this.yieldFunction(stress,q, am, mat);
          if(f <= EP_TOL){
              conv = false;
              if(!conv) throw new ReturnNotAchieved("cutting-plane");
              break;
          }   
          if(k > EP_MAX_ITE){
              conv = false;
              if(!conv) throw new ReturnNotAchieved("cutting-plane");
              break;
          }
          
          IMatrix D = plasticModuli(alpha, mat);
          hm.put(this.dFdSIG,yieldGradStress (stress,q,am,hm));
          hm.put(this.dFdq ,yieldGradIntVar (stress,q,mat,hm));
          v.mul(c,(IVector)hm.get(this.dFdSIG));          
          w.mul(D,(IVector)hm.get(this.dFdq));          
          Double denom = new Double(v.dot((IVector)hm.get(this.dFdSIG)) + w.dot((IVector)hm.get(this.dFdq)));
          hm.put(this.DENOM,denom);
          double ddg = f /((Double)hm.get(this.DENOM)).doubleValue();          
          IVector dep = new IVector(am.getNumberOfDualInternalVariables());
          dep.scale(ddg,(IVector)hm.get(this.dFdSIG));
          ep.add(ep,dep);
          IVector dalpha = new IVector(this.getNumInternalVar(mat));
          dalpha.scale(ddg,(IVector)hm.get(this.dFdq));
          alpha.add(alpha,dalpha);
          Double inc = new Double(((Double)hm.get(this.PLASTIC_MULTIPLIER)).doubleValue()+ ddg);
          hm.put(this.PLASTIC_MULTIPLIER,inc);           
          k++;
        }
        while(true);
        if(conv){
            hm.put(this.ALGORIT_MOD,c);
            this.setPlasticVar(ep, alpha, mat,hm);
        }
        return conv;
    }
        
    /** Stores the plastic strains and the internal variables in the interactive variables.*/
    protected abstract void setPlasticVar(IVector ep, IVector alpha, Material mat, HashMap<Object, Object> hm);
    
    /** Returns the matrix of generalized plastic moduli. 
     *@return The matrix of generalized plastic moduli.
     */
    public abstract IMatrix plasticModuli(IVector alpha, Material mat);
    
    /** Returns the number of internal variables this constitutive model. 
     *@return The number of internal variables.
     */
    public abstract int getNumInternalVar(Material mat);
    
    /** Returns the yield function value from a given state of stress and internal variable. 
     *@return The yield function value.
     */
    public abstract double yieldFunction(IVector s, IVector q, AnalysisModel am, Material mat);    
        
    /** Returns the internal variables in stress space from a given set of variables in strain space. 
     *@return The internal variables in stress space.
     */
    public abstract IVector internalVarStress(IVector alpha, Material mat);
    
    /** Returns the internal variables in strain space this constitutive model. 
     *@return The internal variables in strain space.
     */
    public abstract IVector getInternalVar(Material mat, HashMap hm);
    
    /**Returns the Hessian of yield surface w.r.t. the stress components and internal variables in stress space.
     * @param s The stress vector.
     * @param q The internal variables in stress space. 
     * @return Returns the Hessian of yield surface.
     */
    public abstract IMatrix yieldHessian(IVector s, IVector q, AnalysisModel am);        
        
    /** Returns the gradient of yield surface w.r.t. internal variables in stress space. 
     *@return The gradient of yield surface w.r.t. internal variables in stress space.
     */
    public abstract IVector yieldGradIntVar(IVector s, IVector q, Material mat, HashMap<Object, Object> hm);       
    
    /** Returns the gradient of yield surface w.r.t. the stress components. 
     *@return The gradient of yield surface w.r.t. the stress components.
     */
    public abstract IVector yieldGradStress(IVector s, IVector q, AnalysisModel am, HashMap<Object, Object> hm);
       
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#getNumberOfVariables()
     */
    public abstract int getNumberOfVariables();
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#getNumberOfPreviousVariables()
     */
    public abstract int getNumberOfPreviousVariables();
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#update(java.util.HashMap, java.util.HashMap)
     */
    public abstract void update(HashMap<Object, Object> previous, HashMap<Object, Object> cv);
    
    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.materialMedia.constitutiveModel.ConstitutiveModel#init(java.util.HashMap, java.util.HashMap, br.ufmg.dees.insane.analysisModel.AnalysisModel, br.ufmg.dees.insane.materialMedia.material.Material)
     */
    public abstract void init(HashMap<Object, Object> previous, HashMap<Object, Object> cv, AnalysisModel am, Material mat);
    
    private void check(AnalysisModel am, Material mat){
    	try {
			this.check(am);
			this.check(mat);
		} catch (AnalysisModelNotValid e) {
			e.printStackTrace();
            System.exit(0);
		}catch (MaterialNotValid e) {
			e.printStackTrace();
            System.exit(0);
		}
    }

	public IMatrix mountCt(ContinuousPointModel am, Material mat, HashMap cv) {
		// TODO Auto-generated method stub
		return null;
	}

	public IMatrix mountCs(ContinuousPointModel am, Material mat, HashMap cv) {
		// TODO Auto-generated method stub
		return null;
	}
    
}

final class ReturnNotAchieved extends Exception{
    /**
     * Comment for <code>serialVersionUID</code>
     */
    private static final long serialVersionUID = 1L;
    
    ReturnNotAchieved(String message){
        super("Return not achieved with the" + message + "algorithm");
    }
}
